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We study the evolution of a localized perturbation in a 
chemical system with multiple homogeneous steady states, in 
the presence of stirring by a fluid flow. Two distinct regimes 
are found as the rate of stirring is varied relative to the rate 
of the chemical reaction. When the stirring is fast localized 
perturbations decay towards a spatially homogeneous state. 
When the stirring is slow (or fast reaction) localized pertur- 
bations propagate by advection in form of a filament with a 
roughly constant width and exponentially increasing length. 
The width of the filament depends on the stirring rate and re- 
action rate but is independent of the initial perturbation. We 
investigate this problem numerically in both closed and open 
flow systems and explain the results using a one-dimensional 
"mean-strain" model for the transverse profile of the filament 
that captures the interplay between the propagation of the 
reaction-diffusion front and the stretching due to chaotic ad- 
vection. 



The problem of chemical or biological activity 
in fluid flows has recently become an area of ac- 
tive research Apart from theoretical inter- 
est this problem has a number of industrial |9|] 
and environmental [10,11 1 applications. One of 
the simplest manifestation of non-linear behavior 
in reaction-diffusion systems is the possibility of 
travelling front solutions. In this paper we study 
the effect of chaotic mixing, by an unsteady lami- 
nar flow, in reaction-diffusion systems supporting 
front propagation. 



I. INTRODUCTION 

Let us consider N interacting chemical or bio- 
logical components, with dimensionless concentrations 
Ci(r,t),i = 1, ..,N, advected by an incompressible fluid 
flow, v(r, t), that is assumed to be independent of 
the concentrations. The spatio-temporal dynamics of 
the fields is governed by the set of reaction-advection- 
diffusion equations 
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a + v(r, t)-VC i =F i (C 1 ,..,C N ;k 1 ,.., k M ) 



(1) 



where the functions Ti describe the interactions between 
the different components. These may be chemical re- 
actions or, in the case of biological populations (e.g. 
different plankton species), they may represent growth, 
grazing, reproduction, death, predation etc p2[ , p!3[ . The 
parameters ki are the reaction rates characterizing the 
speed of the chemical or biological interactions and k is 
the diffusivity, assumed to be the same for all compo- 
nents. We assume that the flow is laminar (smooth) and 
time-dependent implying chaotic advection |14|-16|, i.e. 
fluid elements separate exponentially at a rate given by 
the Lyapunov exponent A, of the advection dynam- 
ics. 

The equation (|l|) can be non-dimensionalized by the 
transformations 
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where L and U are the characteristic length-scale and 
velocity of the flow and one of the reaction rates, k = 
ki, is used to define a characteristic chemical time-scale. 
Thus, the non-dimensional problem can be written as 



^C i + v(r,t)-VC i =DcLF i (C 1 ,..,C N ) + Pe- 1 AC i 
at 



where 



kL LU 
Da — — and Pe = 

U K 



(3) 



(4) 



are the Damkohlcr and the Peclet number, respectively. 
The Damkohler number jl^JT^], Da, characterizes the 
ratio between the advective and the chemical (or biologi- 
cal) time-scales. Large Da corresponds to slow stirring or 
equivalently fast chemical reactions and vice versa. The 
Peclet number, Pe, is a measure of the relative strength 
of advective and diffusive transport. We are going to 
consider the case of large Pe number, typical in many 
applications, when advective transport dominates except 
at very small scales. 

In applications it may be useful to understand the be- 
havior of a chemical system for a range of stirring speeds 
when other parameters (e.g. reaction rate constants, dif- 
fusivity) are kept unchanged. Although the stirring speed 
affects both non-dimensional numbers in (Q), its varia- 
tion leaves the family of curves on the Da — Pe plane 
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defined by DaPe = constant invariant. For this case the 
appropriate non-dimensional equation could be obtained 
by dividing both sides of (|J) by Da, implying the use 
of the chemical time-scale, 1/fc, as the time unit. Then 
the two non-dimensional parameters would be Da and 
DaPe = kL 2 / k, where the latter, representing the ratio 
of the diffusive and the chemical time-scales, is indepen- 
dent of the stirring rate. 

In this paper, we study the behavior of the chemical 
system for different Damkohler numbers, when the Peclet 
number is kept fixed. Strictly speaking, this corresponds 
to the variation of the chemical reaction rates and cannot 
be achieved by changing the stirring rate alone (except 
if the diffusivity is also changed in order to keep Pe con- 
stant). 

For simplicity, in the following we consider the case 
of two-dimensional flow only, since it is computationally 
cheaper and easier to visualize. This is directly relevant 
to certain geophysical problems where density stratifica- 
tion makes the flow quasi-two-dimensional (e.g. strato- 
spheric chemistry, plankton in the ocean surface layer) 
and is also accessible to laboratory experiments (using 
soap-films, stratified fluids etc.). However we believe 
that many of the results presented in this paper can be 
straightforwardly extended to three dimensional systems. 

With the above assumptions the problem defined by 
Eq.^ is still very general until the interaction terms, T% 
are specified. Previous work has investigated the spatial 
structure of the chemical fields in the case of chemical dy- 
namics with a single stable local equilibrium concentra- 
tion, that is a function of the spatial coordinate [^0] . 
Here we consider the case of chemical systems that in the 
spatially homogeneous case would have multiple steady 
state solutions. There are many examples of multiple 
steady states in interacting chemical or biological sys- 
tems, in models of atmospheric chemistry ]24| , or in the 
dynamics of plankton populations Jl^,|25| . Perhaps the 
simplest is the auto-catalytic reaction A + B — > 2B. 

Thus we assume that the dynamical system that de- 
scribes the evolution of the spatially uniform chemical 
system, Cj(r,t) = Ci(t), 



dC 

-^=DaF i (C 1 ,..,C N ) 



(5) 



has more than one, stable or unstable, fixed point, or 
equivalently the system of equations 



•FiiCy , .., C* N ) — 



(6) 



has multiple roots in the positive quadrant. (Since the 
concentration fields must be positive everywhere, only 
the fixed points with C* > are relevant.) 

We study the response of the system, initially in one 
of the uniform steady states, to a localized perturbation. 
By localized we mean, that the the spatial extent of the 
perturbation, 5, is much smaller than the characteristic 
length-scale of the velocity field, (!<1. The evolution of 
the chemical fields is investigated for different values of 



the Damkohler number in both closed flow and open flow 
systems. The cases of stable and unstable initial uniform 
states will be considered separately, each represented by 
a simple model system. 

In the following section we describe the models used 
for the reaction dynamics and for the flow, followed by 
the presentation of explicit two-dimensional numerical 
simulations in section III. Then, in section IV. we intro- 
duce and investigate a one-dimensional reduced model 
and show that this may be used successfully to interpret 
the two-dimensional numerical simulations. The paper 
ends with a summary and discussion. 



II. THE MODELS 

For the reaction term, T ', we use two different models 
with multiple equilibria with quadratic and cubic non- 
linearity. The first is an auto-catalytic reaction, a generic 
model for the propagation of a stable phase into an un- 
stable one. The second model is a bi-stable system and 
we study the triggering of a transition from one stable 
state to the other by a localized perturbation. 

1. Auto-catalytic reaction - Let us consider the 
reaction A + B — > 2B, with the corresponding rate equa- 
tions for the spatially uniform system 



dC A 
dt 



-rC A C E 



dC B 
dt 



rC A C E 



(7) 



Observing that C A + Cb (the total number of molecules 
A and B) is conserved by the reaction, we can eliminate 
C A and characterize the state of the system with a new 
variable C = Cb/{C a + Cb) (0 < C < 1) representing 
the proportion of component B, that evolves according 
to 



dC 

^ = kT{C) = kC{l-C), 



(8) 



where k = r(C A + Cb)- 

There are two steady states: i.) C — (component 
A only). This is unstable, since the addition of a small 
amount of B leads to the complete consumption of A 
through the auto-catalytic reaction; ii.) C — 1 (compo- 
nent B only). This is stable - a small amount of A is 
quickly transformed back to B. 

The temporal evolution of the spatially homogeneous 
system can be obtained by integrating m) as 



C(t) 



' 1 C(0) ' 



-i -l 



(9) 



where C(0) is the initial proportion of component B. 

In the numerical experiments, the initial state is the 
unstable phase, C = 0, perturbed by a localized pulse of 
the form 



C(x,y,t = Q) = C e-^+v 2 y^ 



(10) 
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where Co <C 1 and S <C 1, representing a small dilute 
patch of the catalyst, component B. 

The reaction term (|J) can also be interpreted as the, 
so called, logistic growth of a population [2(| modelling 
the growth of the population limited by the availabil- 
ity of the resources (e.g. nutrients), implying that the 
concentration saturates at C = 1. Furthermore, the 
same reaction term appears in the Kolmogorov-Piskunov- 
Petrovsky (KPP) equation |27|] used in combustion P,p8[, 
which describes the propagation of a flame front separat- 
ing fresh (unburned) premixed reactants (C ~ 0) and 
burned gases (C = 1). Thus our numerical experiments 
can also be regarded as representing the evolution of a 
plankton bloom stirred by ocean currents or of a flame 
embedded in a laminar flow. 

2. Bi-stable system - As a second model we use 
a system with two stable states defined by the reaction 
equation 

dC 

— = kF(C) = kC(a- C)(C- 1), < a < 1 . (11) 

The stable fixed points are C = and C = 1, and their 
basins of attraction are separated by the unstable fixed 
point, C = a. Although this system does not represent a 
real chemical reaction scheme, it is the simplest possible 
model for multi-stability. We expect that more realis- 
tic chemical and biological systems with multiple stable 
steady states would show similar behavior. We start the 
system in the stable uniform state, C = 0, and a localized 
perturbation of the form ([Ic]) is added. Now the ampli- 
tude of the perturbation is chosen such to exceed the 
threshold, a, (we use Co = 1.0). (Since the initial state 
is stable, any perturbation that is below the threshold 
everywhere, C(x, 0) < a, dies out.) 

We note that when C = a is chosen as an initial state 
the bi-stable model shows qualitatively similar behavior 
to the auto-catalytic model. 

Stirring will be modelled by two simple time-periodic 
model velocity fields, representing a closed and an open 
flow system, respectively. We stress, however, that the 
results described in this paper are valid for a wide class 
of two-dimensional time-dependent laminar flows, since 
the only important assumption is the chaotic motion of 
fluid elements. 

For the closed flow we choose a sinusoidal shear flow 
with the direction of the shear periodically alternating 
along the x and y axis |^9|,|o) . The velocity field is 

v x (x,y,t) = AQ(^ - t mod l) sin (2i:y + fa) 

v y (x, y, t) = A<d(t mod 1 — sin (2nx + fa+i), (12) 

defined on a doubly periodic square domain of unit 
length, O is the Heaviside step function and A is a pa- 
rameter (we use A = 0.7), that controls the chaotic be- 
havior of the flow. To avoid transport barriers (due to 
KAM tori |0, typically present in periodically driven 



conservative systems) the periodicity is broken by using 
a random phase, fa, different in each time period. 

We also consider open flow systems in which fluid con- 
tinuously flows in and out a finite mixing zone. Such sys- 
tems are relevant for certain chemical reactors and also 
in some geophysical problems. Advection in this type of 
open flows has been shown to be governed by a chaotic 
scattering type escape process generating fractal patterns 
of the advected tracers [ pl|j3^ ] . 

As an example of an open flow system we use a ve- 
locity field modelling the flow around two alternately 
opened point-sinks in an unbounded two-dimensional do- 
main The fluid particles approach the mixing 
zone from infinity and leave the domain through one of 
the sinks. The velocity field is composed by the superpo- 
sition of a point-vortex and a point-sink, centred on the 
active sink. The complex potential corresponding to the 
vortex-sink is 



,(z) = -(Q + iK)ln\z-z s \ 



(13) 



where z = x + iy is the complex coordinate, z s is the po- 
sition of the sink and Q and K are the sink-strength and 
vortex-strength, respectively. The corresponding veloc- 
ity field in polar coordinates with the origin fixed to the 
active sink is 



Q 



K 



v r = — -, , V4, = — , r = y/(x- x s ) 2 + (y- y s ) 2 - 
r r 

(14) 

The half distance separating the two sinks is assumed 
to be unity (x s = 0,y s — ±1). The sinks are alter- 
nately opened for equal times and the period of the flow 
is the time unit. The inflow concentration at the bound- 
aries of a square domain is kept constant at the value 
corresponding to the initial background concentration, 
C(±3.0,y) = C(x,±3.0) = 0.0. 

For both closed and open flows we integrate 
the advection-reaction-diffusion problem for the auto- 
catalytic and the bi-stable model on a 1000 x 1000 lattice 
using a semi-Lagrangian scheme. The value of the Peclet 
number is Pe = 10, 000, and the Damkohler number is 
varied in a range from zero to few hundred. 



III. NUMERICAL RESULTS 
A. Closed flow 

We find for both chemical models two distinct regimes 
separated by a critical Damkohler number, Da c . The 
critical values are Da c ~ 2.0 for the auto-catalytic reac- 
tion and Da c w 20.0 for the bi-stable model. 

In the slow reaction/fast stirring regime (Da < Da c ), 
in both models the initial perturbation quickly decays 
towards a homogeneous state (Fig.l). In the case of the 



3 



auto-catalytic reaction the homogenisation of the per- 
turbation is followed by a spatially uniform transition 
to the stable equilibrium, C = 1, as in a reactor with 
initially premixed components. This is because the orig- 
inal uniform state is unstable, and cannot be restored 
by the homogenisation of the perturbation, since the ho- 
mogenised state still deviates slightly from the unstable 
equilibrium. In the bi-stable system the perturbation is 
dispersed and the system simply returns to the unper- 
turbed initial state, C = 0. In this regime, the chemical 
reaction is too slow to sustain the localized perturbation 
that is diluted by the strong stirring. Note that for both 
chemical models the final state would be the same for 
a spatially uniform perturbation with the same spatial 
mean. Thus, a coarse grained model could, at least qual- 
itatively, reproduce the evolution of the system. 

When Da > Da c the localized perturbation may per- 
sist and propagate in the form of a filament with a 
roughly constant width and rapidly increasing length 
(Fig. 2). This continues until the filaments cover the 
whole domain and finally the system becomes uniform 
again, C = 1. This occurs in all cases in the auto- 
catalytic case, but only if a < 0.5 in the bistable case. If 
a > 0.5 a localised perturbation in the bistable case de- 
cays. The width of the propagating filaments increases 
with Da. We note, that the average profile of the fil- 
ament (width, amplitude) is apparently independent of 
the details of the initial perturbation, indicating that it 
is determined by the interplay between the chemical and 
transport dynamics in the system. In this regime the 
transition from the initially uniform state to the final 
one is strongly non-uniform in space. Since the spatial 
variation is essential, a coarse grained model that was 
unable to resolve the filaments, would produce very dif- 
ferent outcomes. 

The transition from the non-homogeneous to homoge- 
neous reaction in the case of the auto-catalytic model 
for decreasing Da is clearly shown by the snapshots of 
the spatial distribution taken at the midpoint of the 
transition defined by the spatial mean concentration 
C = 0.5 (Fig. 3.). To characterize the change in the non- 
uniformity of the reactions as Da is varied we plotted the 
relationship between the first and the second moments 
(Mi = (C) and M 2 = (C 2 )) of the chemical distribu- 
tion (Fig. 4). There are two extreme situations. In case 
of a spatially uniform system the averaging can be ig- 
nored and thus M 2 = M\. For a strongly non- uniform 
distribution with only two possible values 1 and (repre- 
senting a two-phase system with a very narrow transition 
zone) the square for the second moment is irrelevant and 
Mi = Mi. Fig. 4. clearly shows the transition from the 
linear to the quadratic relationship as Da is decreased. 

The time-dependence of the mean concentration, C, 
(equivalent here to the spatial average of the deviation 
from the initial state) is shown in Fig. 5 for different Da 
numbers. For large Da numbers, after an initial transient 
time, an exponential growth can be observed with the 
growth rate independent of the Damkohler number for 



both chemical models. This shows that the growth of the 
mean concentration is controlled by the stirring, that in- 
creases the length of the propagating filament. Note, that 
for the bi-stable system a localized perturbation with its 
spatial mean concentration well below the threshold, a, 
can flip the system to the other steady state. This is a 
strong example where spatial smoothing does not work. 

For Da < Da c the homogeneous dynamics is relevant. 
In the bi-stable system the mean concentration simply 
decays to zero exponentially as in the homogeneous sys- 
tem (for C < 1) 



dC dT- 

H = Da dc c 



-uDaC, C 



-aDat 



(15) 



In the auto-catalytic model, the growth of the mean con- 
centration depends on the Da number. When the time- 
dependence is plotted against Dot (Fig. 6), (i.e. the time 
unit is the chemical time) the curves corresponding to 
Da < Da c collapse showing that in this regime the tran- 
sition is independent of the stirring rate, as expected for a 
spatially uniform system, and the numerical results agree 
well with the solution obtained for the homogeneous sys- 
tem (^). In this regime the two reactants (A and B) 
are brought close to each other by the flow at a higher 
rate than they can react, therefore further increase of the 
mixing rate cannot enhance the production. The growth 
of the mean concentration is limited by mixing (trans- 
port) for supercritical Damkohler numbers and limited 
by the chemical reaction (local dynamics) below the crit- 
ical value. 



B. Open flow 

In the open flow case we find again a transition at a 
critical value of the Damkohler number. 

When the stirring is strong (Da < Da c ) any local- 
ized perturbation dies out and both models return to the 
original state. There is no homogeneous transition for 
the auto-catalytic system as this would be incompatible 
with the inflow boundary conditions. In other words, 
the perturbation is expelled completely from the system 
(through the sinks) and thus even the unstable basic state 
(C = 0) can be restored. 

When Da > Da c , similarly to the closed flow case, 
the perturbation produces a propagating filament. How- 
ever, in this case due to the continuous outflow and in- 
flow the filament cannot fill the domain uniformly. Af- 
ter a short transient a non-uniform stationary state sets 
in, with the mixing zone partly covered by a complex 
filamental structure (Fig. 7). In our case the stationary 
pattern varies periodically with the period of the flow. 
For the auto-catalytic model, a small amplitude pertur- 
bation is sufficient to initiate the propagating filament, 
while in the bi-stable model only perturbations larger 
than the threshold, a, are able to trigger the transition 
to the non-uniform stationary state. (Furthermore an 



4 



additional condition for the existence of the non-uniform 
stationary state is again that a < 0.5.) 

One can observe that the non-uniform (periodic) sta- 
tionary pattern follows the fractal unstable manifold of 
the non-escaping set formed by fluid particles that are 
trapped forever in the mixing zone |34| . The unstable 
manifold can be easily visualized by simply following the 
evolution of an ensemble of fluid elements (e.g. represent- 
ing a droplet of dye) advected by the flow (Fig. 8). For the 
auto-catalytic reaction this has been already observed in 
a kinematic model where B particles are treated as indi- 
vidual tracers |55| . As in the closed flow case, the width 
of the filaments increases with the Da number (Fig. 9). 
The dependence on the Da number of the total concen- 
tration (i.e. the area covered by the filaments) in the 
stationary state is shown in Fig. 10. We find a continuous 
transition for the auto-catalytic reaction (Da c = 2.3) and 
a discontinuous one for the bi-stable model (Da c = 24.2). 

IV. THE LAGRANGIAN FILAMENT SLICE 
MODEL 

Here we introduce a reduced one-dimensional model 
in order to explain the numerical results described in 
the previous section. In the presence of chaotic trans- 
port fluid elements are stretched into elongated fila- 
ments. This is well known from numerical simulations 
and has been observed in laboratory experiments using 
dye droplets p6| , ^2|j3^ ] . In a two-dimensional system, one 
can assign to any point of the flow a convergent and a 
divergent direction associated with the eigenvectors cor- 
responding to the negative and positive Lyapunov expo- 
nents (—A and A) of the chaotic advection. These direc- 
tions are tangent to the stable and unstable foliations of 
the advection dynamics, respectively. Any advected ma- 
terial line (e.g. iso-contours of a conserved tracer) tends 
to align along the unstable foliation in forward time, or 
along the stable foliation in backward time. Thus, the 
stirring process smooths out the concentration of the ad- 
vected tracer along the stretching direction, whilst en- 
hancing the concentration gradients in the convergent 
direction. 

Let us now separate the original reaction-advection- 
diffusion problem along the (Lagrangian) stretching and 
convergent directions. In the stretching direction the 
perturbation is spread by advective transport, that is 
the dominant process being much faster than diffusion. 
In the convergent direction the formation of small scale 
structure indicates that all three processes of reaction, 
advection and diffusion are important and need to be 
considered together. The resulting equation determines 
the mean transverse profile of the filament that propa- 
gates along the divergent direction following the unstable 
foliation. 

Thus one expects that the locus of the centre of the 
filament can be obtained simply by advecting a material 



contour starting in the region of the initial perturbation. 
This is confirmed by the numerical simulations as shown 
in Fig. 11 for the closed flow model, and is consistent with 
Fig. 8 for the open flow system. An important difference 
between the two is that while in the closed flow the con- 
tour gradually fills the whole domain, in the open flow it 
draws out the unstable manifold of the set formed by all 
non-escaping orbits. The length of the contour increases 
exponentially (Fig.12), L(t) ~ exp (9t) with 9 w 2.05. 
We note, that the contour lengthening rate, 9 is always 
larger than the Lyapunov exponent (A), which represents 
the average growth rate of infinitesimal line elements. 
This is because the instantaneous stretching rate fluctu- 
ates and the increase of the total length is dominated 
by the growth of a line elements that experience a faster 
than average stretching. In dynamical systems language 
the contour lengthening rate 9 is given by the topological 
entropy |l7||37|,[$8| of the advection dynamics. 

In the convergent direction we have the following one- 
dimensional equation for the average profile of the fila- 
ment 

® C -\ x ^-C = kT{C) + K^C, (16) 
at ox dx 

representing the evolution of a transverse slice of the fil- 
ament in a Lagrangian reference frame (i.e. following the 
motion of a fluid element). The second term on the left 
hand side is a stretching term that takes into account the 
mean convergent flow and T is the original reaction term. 
The stretching term in equation ( [l6| ) can be interpreted 
as advection by a pure strain flow along its convergent 
direction v x = —\x. In the two-dimensional problem the 
strength of stretching fluctuates in space and time. To 
capture the average behavior this can be represented by 
a constant stretching rate, A, equal to the Lyapunov ex- 
ponent of the chaotic advection. Therefore the equation 
above can be regarded as a Lagrangian mean field de- 
scription. Equation (jl^) has also been studied recently 
by McLeod et al investigating the filament width of 
oceanic plankton distributions. 

The equation ([l6]) is defined on — oo < x < oo with 
the boundary conditions 

dC 

C(x -> ±oo) = 0, — (a: ±oo) = 0. (17) 
ax 

representing the assumption that most of the system is 
in the background state, C = 0, so that different parts 
of the filament are well separated from each other and 
they do not interact. The single filament approximation 
is not valid for the late stage of the evolution when the 
filaments can overlap. 

The homogeneous steady state, C(x) = 0, is a trivial 
solution of ( ^6[]l7| ). Let us now consider the evolution of 
a localized pulse-like disturbance. (We use perturbations 
centred at the origin, but it turns out that the asymp- 
totic behavior is independent of the initial position of the 
perturbation.) 
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For a non-reactive tracer (k = 0) the equation Jl6| ) 
has asymptotic solutions (for large t) of the form of a 
Gaussian pulse whose amplitude decays exponentially in 
time 



C(x, t) ~ exp {—Xt) exp 



x 2 X 
~2k~ 



(18) 



The width of the Gaussian, Idif = \J n/X, is determined 
by the balance between the strain and diffusion. 

Let us now consider the case of a reactive tracer 
(k 0) without stretching (A = 0). It is well known 
that reaction-diffusion systems with multiple equilibria 
have travelling front solutions connecting different steady 
states |^6|,[40],[n] . The fronts move with a constant speed 
vq with no change of shape, C{x,t) = C(x — vot). 

The reaction-diffusion problem corresponding to the 
auto-catalytic model is known as the Fisher equation |^] 
(or Kolmogorov-Petrovsky-Piskunov equation |27j in the 
combustion literature) and describes the propagation of a 
stable phase (C = 1, component B) into an unstable one 
(C = 0, component A). Localized perturbations evolve 
into a pair of fronts moving away from the centre with 
the asymptotic speed, vo = 2\fkn. 

For the bi-stable model ([ll]) the velocity of the front 
joining the two stable states, C(x — > — oo) = 1 and 
C(x -> oo) = 0, is pGlMplJ 



t'o 



HskJ^ T{C)dC = (1 ~ 2a). 



(19) 



The sign of the above expression can be either positive 
or negative showing that the direction of the propagation 
depends on the parameter a. The single travelling front 
solution can be found analytically pq] as 



C(x — vot) 



1 + exp 



x — Vot 

~7T 



(20) 



When the initial basic state is C — a localized pertur- 
bation can initiate a pair of propagating fronts moving 
away from each other only if a < 0.5, otherwise the di- 
rection of the front propagation is such that the fronts 
approach each other and the perturbation dies out. This 
explains the decay of the perturbations for a > 0.5 in the 
two-dimensional simulations independently of the stir- 
ring rate. 

Thus, in the absence of stretching both type of sys- 
tems have travelling front solutions with the front speed 
proportional to \fk~K,. A localized perturbation generates 
a pair of fronts moving in opposite directions away from 
the centre. (For bi-stable systems this happens only if 
a < 0.5). For the auto-catalytic model the amplitude of 
the perturbation can be arbitrarily small, while it must 
exceed the threshold a in the bi-stable case. 

In the presence of stretching, A > one expects that 
the convergent flow will eventually stop the propagation 
of the fronts at a point x = w where the speed of the 
propagation is balanced by the advection 



wX 



'kn 



(21) 



This gives an estimate for the width of the resulting fil- 
ament as 



Da = 



A 



(22) 



where we have introduced the Lagrangian Damkohler 
number, Da. This is defined on the basis of the Lya- 
punov time, 1/A, of the flow, representing the Lagrangian 
characteristic time-scale of the two-dimensional advec- 
tion problem. This turns out to be more appropriate for 
the filamentation problem than the definition (Q) based 
on Eulerian characteristics like the average flow velocity 
U. 

The propagation velocity (]2l]) is for an isolated front 
only. Therefore we expect that the scaling for the fil- 
ament width ( p2| ) is valid when the distance between 
the two fronts, representing the edges of the steady fila- 
ment, is sufficiently large compared to the diffusive scale, 
w S> Idif, that is Da ^> 1. 

Equation ( |l6| ) can be non-dimensionalized by using the 
Lyapunov time, A" 1 , as the time-scale unit and the dif- 
fusive scale, Idif, as the unit length 



c) 



d 



d 2 



^C-x-C = Da^C) + -^C. 



dt 



dx 



dx z 



(23) 



Since the one-dimensional problem is defined on an un- 
bounded domain, the Peclet number does not appear in 
(p3|). For the two-dimensional problem the characteris- 
tic scale of the velocity field, L, is finite and this can be 
used to define a Lagrangian Peclet number based on the 
Lyapunov exponent of the flow as 



Pe = 



L 2 X 



Idif J 



(24) 



Using this the expression for the width of the filament 
can be rewritten as 



L 



XL 




(25) 



The straining flow approximation can only be valid for 
scales smaller than L (w <C L), thus (^5|) is expected to 
be correct in the range 1 <§; Da <§; Pe. 

In the open flow model in the stationary state the fila- 
ments cover a fractal set. For the parameter values used 
in our simulations the dimension of this was found to be 
D w 1.69. The area of the large concentration region can 
be obtained by using the dimensionless filament width 
w/L as the resolution when observing this region. Since 
the number of boxes of size w/L needed to cover the re- 
gion is proportional to (w/L)~ D |rij] , the area A of this 
regions is 



A 
I* 



U)\2- Z> 

~L, 




(2-D)/2 



(26) 



G 



This shows, that due to the overlaps the total area grows 
more slowly with the Damkohler number than the area of 
a single isolated filament. An analogous scaling has been 
obtained for auto-catalytic reaction on an open baker 
map in |43|. 

Numerical simulations of the reduced problem (|23| ) 
show that for both chemical models there is a critical 
value of the Damkohler number. When Da > Da c there 
exists a non-uniform steady solution to ( |23| ) centred on 
the origin (Fig. 13). Otherwise all perturbations decay 
and the only steady solution is the trivial one, C(x) = 0. 
The width of the non-uniform steady solution increases 
with Da, that appears to be consistent with (|2^) for large 
Damkohler numbers. Numerical continuation of the non- 
uniform steady solution for decreasing Da confirms that 
this solution disappears at a critical value. The transition 
is found to be continuous for the auto-catalytic model and 
discontinuous for the bi-stable model (see Fig. 13). Let us 
now analyse the transition in the reduced problem ( [23] ) 
in more details for the two chemical models, separately. 



A. Auto-catalytic model 

For the auto-catalytic model the non-uniform solution 
approaches and coalesces with the uniform solution when 
Da c is approached from above. When Da < Da c local- 
ized perturbations decay, showing that the uniform state, 
C(x) = 0, in spite of being an unstable fixed point of 
the homogeneous system is stable against localized dis- 
turbances. Thus stirring stabilizes the unstable equilib- 
rium of the homogeneous system by dispersing and di- 
luting the perturbation. This is consistent with the be- 
havior observed in the two-dimensional simulations for 
sub-critical Damkohler numbers, showing homogenisa- 
tion and decay of the perturbation followed by growth 
only after the reactants were distributed uniformly in 
space. For supercritical Damkohler numbers the uniform 
solution is unstable against localized perturbations. Ar- 
bitrarily weak perturbations are sufficient to reach the 
non-uniform steady state. 

Just above the critical point, Da—Da c <C 1, the ampli- 
tude of the non-uniform solution is small and the chem- 
ical dynamics can be linearized about the background 
state 



dC 
~dt 



dC 
dx 



DaC 



d 2 C 
dx 2 



This problem has been studied by Martin 
text of plankton populations. Equation (E 
totic solutions of the form 



C{x,t) 



1)' 



(27) 

I in the con- 
has asymp- 

(28) 



that decay in time when Da < 1 and grows exponentially, 
moving out of the domain of validity of the linear approx- 
imation, otherwise. (The non-linearity would eventually 



stop the growth.) Thus we obtain the critical value for 
the auto-catalytic model, Da c = 1. 

Close to the transition, Da = l + e,0<e<Sl ) one can 
look for a steady non-uniform solution of the form 



C{x;e) = eC\{x) + e 2 C 2 (x) + ... 



(29) 



Substituting this into ( jig ) for the term first order in e we 
obtain 



d 2 C, dCi 

-TT + x ~j ^1=0 

ax* ax 



(30) 



that has the solution C\ = Ae~ x / 2 , where A is a con- 
stant that can be determined from the equation for the 
terms second order in e 



d 2 C 2 , dC 2 



C 2 = C^-d. 



(31) 



dx 2 dx 

The left hand side of the equation can be written as 

Integrating both sides over the whole domain the left 
hand side vanishes, according to (|l7|), and an equation 
for the constant A is obtained 



I" (C 2 - d)dx = A 2 ^ - A^- 



(33) 



that gives A = l/\/2. Thus, the steady solution close to 
the transition point can be approximated as 



C(x;e) 



v/2 



0{e 2 ) 



(34) 



B. Bi-stable model 

For the bistable case the transition at Da c is discon- 
tinuous. The non-uniform solution disappears with a fi- 
nite amplitude far from the uniform state. The uniform 
solution is stable for any Da thus small perturbations 
decay independently of the Damkohler number. In the 
supercritical regime the uniform and non-uniform stable 
solutions coexist suggesting the presence of a threshold 
for exciting the non-uniform perturbation. 

To investigate the transition further let us look for 
steady solutions of eq. (p3|) 



(35) 



d 2 C ~ dC 
— = -DaT{C) - x— 
ax z ax 



that are consistent with the boundary conditions (Jl_7 
Since the non-uniform steady state is symmetric about 
the origin it is sufficient to consider the domain < x < 
oo, with the constraints 



7 



C(x->oo) -*0,— (0) = 0. (36) 
ax 

Note that, if x is interpreted as time and C as a spa- 
tial coordinate, the problem ( p5| ) is equivalent to the mo- 
tion of a particle in an asymmetric two-humped potential 
(Fig. 14) 

V(C) = -£aC 2 Qc 2 + { -±±^C - |) , (37) 

under the effect of linear friction, with friction coefficient 
increasing linearly in time. The two maxima of the po- 
tential are at the stable fixed points C — and C = 1 
and the minima is at C — a. The potential difference be- 
tween the two maxima is proportional to Da. The parti- 
cle trajectory satisfying (|3(]) starts from the left slope of 
the higher potential 'hill (a < C(x = 0) < 1) with zero 
velocity and ends exactly on the top of the lower hill. 
Thus the problem reduces to finding the appropriate val- 
ues of the initial coordinate, Cq = C(x = 0). For initial 
conditions Cq G (a,l),C'(x = 0) = the trajectory of 
the particle may either end in the potential well C = a 
or may cross the smaller hill and escape to — oo. The 
trajectories corresponding to the non-uniform steady so- 
lutions are at the boundary between these two types of 
asymptotic behavior. 

We calculated numerically trajectories for a set of ini- 
tial conditions in the range Co G (a, 1) for a set of dif- 
ferent values of the Damkohler number, Da <G (0,40). 
The asymptotic behavior of the trajectories is indicated 
on the Co — Da plane (Fig. 15): blank, C(x — > oo) — > a; 
black, C(x — > oo) — > — oo. The required steady solutions 
are on the boundary of the two regions. The numerical 
results show that for small Da there is no such bound- 
ary and a solution of the type ( |36| ) does not exist. In 
the particle analogy the interpretation of this is that the 
difference in the height of the two hills is not sufficiently 
large to compensate for the energy dissipated by friction, 
thus the particles cannot escape. When Da is increased 
above the critical value Da c there are two solutions. If 
the hill at C = 1 is high enough there are initial con- 
ditions for which the particles have sufficient energy to 
cross the potential barrier. Clearly, particles with ini- 
tial conditions below the height of the lower potential 
hill are still trapped in the potential well. Also, particles 
started from a point very close to the top of the higher 
hill are unable to escape since they may spend very long 
time in the neighbourhood of the stationary point and go 
through the potential well at a late time when the fric- 
tion is strong. Thus, for Da < Da c the initial conditions 
^escape ^ wn j cn particles escape to infinity are in an 
the interval of the form 

a < Cas(Da) < C escape < Cv. 2 {Da) < 1 (38) 

where Co,i(-Da) are Co, 2 (-Da) are initial conditions cor- 
responding to steady non- uniform solutions of (|l6|). As 



the Damkohler number is decreased the two solutions ap- 
proach each other and disappear at Da c = 11.0. 

The trajectory starting from C(x = 0) = Co, 2 (-Da) 
clearly corresponds to the steady non-uniform solu- 
tion found in the numerical simulations of the time- 
dependent problem. The solution corresponding to the 
lower branch, C(x = 0) = Co.i(-Da), however, is not 
found as an attractor of the time-dependent problem. 
This suggests that this is an unstable solution playing 
the role of a threshold separating the basins of attrac- 
tion of the uniform and non-uniform stable solutions. (It 
can be shown that all initial conditions that are above 
(below) this separating solution everywhere, converge to 
the non-uniform (uniform) steady state. This, of course, 
does not say anything about initial conditions partly be- 
low and partly above the separating solution.) 

V. SUMMARY AND DISCUSSION 

The one-dimensional Lagrangian filament model jl^ ) 
clearly explains the qualitative features of the two- 
dimensional numerical results. It shows how a steady 
filament profile can arise as a result of the interaction 
between the propagation of a reaction-diffusion front and 
stretching due to chaotic advection. The disappearance 
of the filament type solution for sub-critical Damkohler 
numbers explains the transition observed in the two- 
dimensional simulations. The advective propagation of 
the filament along the unstable foliation of the chaotic 
advection explains the exponential growth of the mean 
concentration. In principle, these results could be used 
as a numerical technique for obtaining the approximate 
spatial distribution of the chemical components by com- 
bining a two-dimensional contour-advection calculation 
with the information about the filament width obtained 
from the steady solution of the one-dimensional model. 

In order to compare the critical Damkohler numbers 
predicted by the one-dimensional mean strain model with 
the ones obtained from the direct numerical simulations 
we calculate the Lagrangian Damkohler numbers cor- 
responding to the two-dimensional problem. The Lya- 
punov exponents of the advection in the two model flows 
was found to be: X c iosed ~ 1-66 and X op en ~ 2.19. 
The critical values of the Lagrangian Damkohler num- 
ber based on these Lyapunov exponents are presented in 
Table 1. and show a very good agreement. 

In our analysis we neglected the fluctuations of the 
stretching rate. In reality there is a distribution of 
stretching rates. The effect of this is visible in the numer- 
ical simulations showing that the width of the filament 
slightly fluctuates in space and time. Also the direction 
of the stretching fluctuates and foldings of the filament 
may lead to large curvatures whose effect is not captured 
by our one-dimensional description. Another effect is the 
non-uniform density of the unstable foliation pointed out 
by Alvarez et al p8[ |. Thus the advected filament can 
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overlap with itself well before it fills the whole domain. 
Some regions of the flow are filled while others are still 
empty. 

Here we investigated only reactive systems described 
by the distribution of a single species. We expect, how- 
ever, that the basic phenomena described in this paper 
remains valid for multi-component reactive systems that 
may have a number of different chemical time-scales. One 
example of this type is the case of excitable systems 
1 26, pl|j45| under stirring by a chaotic flow discussed in 
[ 46 1 . Excitable systems have two different time-scales - 
corresponding to fast and slow components. Although 
these systems only have a single (stable) steady state, 
the 'rest state', they also have a 'meta-stable' excited 
state that persists for a finite time only. Excitable sys- 
tems under stirring exhibit similar behavior to the one 
presented here, including advective propagation in form 
of a steady excited filament and the existence of a crit- 
ical Damkohler number. The one-dimensional excitable 
reaction-diffusion systems have travelling pulse solutions, 
that in the presence of stretching leads to the appearance 
of a steady excited filament solution. This can be simple 
unimodal, as in our case, but it can also have a bi-modal 
structure with the central part of the filament returned 
to the rest state. Thus the existence of an extra chem- 
ical time-scale in this system allows for somewhat more 
complex structures and a further transition in the large 
Da number range, being a transition from coherent to 
non-coherent excitation of the system. 

A nice example of an advectively propagating pertur- 
bation, of the kind described in this paper, has been ob- 
served recently in a so-called 'ocean fertilization' experi- 
ment , where the injection of a trace element affecting 
the plankton ecosystem triggered a phytoplankton bloom 
in the form of an elongated filament, observed on satel- 
lite images. The response of plankton ecosystem models 
to such perturbation in the presence of stirring has been 
studied in . We suggest that similar phenomena could 
also be investigated in laboratory experiments. 
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FIG. 1. Snapshots of the spatial distribution for the 
autocatalytic reaction for Da = 1.0, (< Da c ) at 
t = 0.0, 2.0, 4.0, 6.0, 8.0. The amplitude of the initial pertur- 
bation is chosen to be Co = 0.5 in order to make the initial 
decay visible. 



FIG. 2. Snapshots of the spatial distribution for the au- 
tocatalytic model for a supercritical Damkohler number, 
Da = 7.0 at t = 0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5. The am- 
plitude of the initial perturbation is Co = 0.5. The bistable 
model shows qualitatively similar behavior. 



FIG. 3. Snapshots of the spatial distribution for the auto- 
catalytic model at the midpoint of the transition (defined by 
< C > (t) = 0.5) for Da = 35.0, 12.0, 8.40, 4.1, 2.9, 1.0. 



FIG. 4. The relationship between the first and the second 
moment of the spatial distribution for different values of Da. 
(a) autocatalytic, (b) bi-stable. 
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TABLE I. The critical values of the Lagrangian 
Damkohler number for the two-dimensional simulations and 
for the one- dimensional single filament model 



FIG. 5. The time dependence of the total concentration 
for different values of Da, (a) autocatalytic reaction and (b) 
bistable model. The dashed line indicates the rate of growth 
of the length of a filament due to advection. 

FIG. 6. Same es Fig5. with rescaled time. The dashed 
line shows the time-dependence for the spatially homogeneous 
system. 



FIG. 7. Snapshots of the spatial distribution for the 
autocatalytic reaction in the open flow system for a su- 
percritical Damkohler number, Da = 14.0 > Da c at 
t = 0.0,1.0,2.0,3.0,4.0,5.0. The distribution at t = 5.0 
has already reached the time-periodic stationary state. The 
amplitude of the perturbation is Co = 0.5. For subcritical 
Damkohler numbers the perturbation dies out. (The bi-stable 
model shows similar behavior but for different values of Da.) 



FIG. 8. The evolution of an ensemble of particles (e.g. 
representing a droplet of dye) in the open flow system 
(t = 0.0,2.0,4.0,6.0). 



FIG. 9. Spatial distribution in the open flow for the au- 
tocatalytic model in the time-periodic stationary state for 
Da = 70.0, 24.0,11.8,4.0. 
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FIG. 10. The total concentration in the steady state for 
the open flow system as a function of the Damkohler number 
for the autocatalytic reaction (a) and the bi-stable model (b). 
The dashed line in the inset indicates the predicted asymp- 
totic behavior C to tai ~ Da' 2_D " 2 where D is the fractal 
dimension of the unstable foliation of the non-escaping set, 
D = 1.69 

FIG. 11. Temporal evolution of a material line ad- 
vected by the closed flow. The radius of the initial cir- 
cular contour is r = 0.06 and the figures correspond to 
t — 0.0, 0.5, 1.0, 1.5, 2.0, 2.5 as in Fig. 2, for comparison. 

FIG. 12. The growth of the length of the contour shown in 
Fig. 11 as a function of time. The continuous line represents 
an exponential growth ~ exp (2.05i). 

FIG. 13. Steady solutions of the one-dimensional reac- 
tion- ad vection-diffusion problem, for different values of the 
Damkohler number (a. Da = 80,40,20,10,5.0,2.5,1.25; c. 
Da = 800, 400, 200, 100, 50, 25, 12.5, 11.0) and the dependence 
of the total concentration in the steady state as a func- 
tion of Da, (a,b) autocatalytic, (c,d) bi-stable. The dashed 
line in the inset indicates the predicted asymptotic behavior 

Ctotai ~ V 1 Da. 

FIG. 14. The two-humped potential V(C) for Da = 1 and 
a = 0.25. 

FIG. 15. The shaded area shows initial conditions result- 
ing in the escape of the particle from the potential well. The 
boundary of the shaded area corresponds to the steady fila- 
ments solutions. Note, that the filament solution disappears 
for sub-critical Damkohler numbers (Da = 11.0) 
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